Trapping of Dust by Coherent Vortices in the Solar Nebula 
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Abstract 



We develop the idea proposed by Barge & Sommeria (1995) and Tanga et al. (1996) 
that large-scale vortices present in the solar nebula can concentrate dust particles and 
' facilitate the formation of planetesimals and planets. We introduce an exact vortex 

solution of the incompressible 2D Euler equation and study the motion of dust particles 
in that vortex. In particular, we derive analytical expressions for the capture time and 
the mass capture rate as a function of the friction parameter. Then, we study how 
small-scale turbulent fluctuations affect the motion of the particles in the vortex and 
determine their rate of escape by solving a problem of quantum mechanics. We apply 
these results to the solar nebula and find that the capture is optimum near Jupiter's 
orbit (as noticed already by Barge & Sommeria 1995) but also in the Earth region. 
This second optimum corresponds to the transition between the Epstein and the Stokes 
regime which takes place, for relevant particles, at the separation between telluric and 
$_i ' giant planets (i.e near the asteroid belt) . At these locations, the particles are efficiently 

captured and concentrated by the vortices and can undergo gravitational collapse to 
form the planetesimals. 



1 Introduction 



Many astrophysical objects, ranging from young stars to massive black holes, are surrounded 
by widespread gaseous disks. The existence of a primordial disk around the sun was conjec- 
tured by Kant (1755) and Laplace (1796) in the lS th century to explain the quasi-circular, 
coplanar and prograd motion of the planets in the solar system. Such protoplanetary disks 
have recently been observed with the Hubble Space Telescope in the Orion nebula around 
stars less than one million years old. These gaseous disks can be considered as a by-product 
of the star formation: after the gravitational collapse of a molecular cloud, about 99% of 
the initial angular momentum is spread in an extended disk, while 99% of the mass forms 
the star, whose internal structure is hardly affected by rotation. 

Whenever it has been possible to observe rotating, turbulent fluids with good resolu- 
tion, it has been seen that individual, intense vortices form (Bengston &: Lighthill 1982; 
Hopfinger et al. 1983; Dowling and Spiegel 1990). One of the most striking example is 
Jupiter's Great Red Spot, a huge vortex persisting for more than three centuries in the 
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upper atmosphere of the planet. These coherent vortices are well reproduced in numerical 
simulations (McWilliams 1990, Marcus 1990) and laboratory experiments (Van Heist & Flor 
1989, Sommeria, Meyers & Swinney 1988) of two-dimensional turbulence and their organi- 
zation can be explained in terms of statistical mechanics (Miller 1990, Robert & Sommeria 

1991, Sommeria et al. 1991, Chavanis & Sommeria 1998) Q. It seems therefore natural to 
expect their presence in accretion disks also (Dowling & Spiegel 1990, Abramowicz et al. 

1992, Adams & Watkins 1995). 

However, accretion disks are exceptional among rotating turbulent objects in the strong 
shears that these bodies are believed to possess and this shear might lead to rapid destruc- 
tion of any structures that tend to form. This objection has been overruled by the numer- 
ical simulations of a two-dimensional flow in an external Keplerian shear by Bracco et al. 
(1998,1999) for an incompressible flow and Godon & Livio (1999a, b,c) for a compressible 
flow. Although cyclonic fluctuations are rapidely elongated and destroyed by the shear, 
anticyclonic vortices form and persist for a long time before being ultimately dissipated by 
viscosity. Naturally, this does not prove that coherent structures must form on disks, but 
this strengthens the argument that disks are likely to follow the norm of rotating, turbulent 
bodies. Other numerical results (Hunter & Horak 1983) and experimental work (Nezlin & 
Snezhkin 1993) comfort this point. 

Coherent vortices in circumstellar disks can play an important role in the transport of 
dust particles and in the process of planet formation. Planets are thought to be formed 
from the dust grains embedded in the disk after a three-stage process: (i) in a first stage, 
microscopic particles suspended in the gas stick together on contact due to electrostatic or 
surface forces. When they reach sufficient sizes, they begin to sediment in the mid-plane 
of the disk due to the combined effect of the gravity and the friction with the gas. When 
settling dominates, a particle can grow by sweeping up smaller ones (Safronov 1969) and 
may easily reach sizes of several centimeters in a few thousand orbital periods (Weiden- 
schilling & Cuzzi 1993). Bigger aggregates (> 100cm) are more difficult to form on relevant 
time scales because of collisional fragmentation (ii) When the layer of sedimented particles 
is sufficiently dense, the gravitational instability is triggered and the layer crumbles into 
numerous kilometer-sized bodies, the so-called "planetesimals" (Safronov 1969, Goldreich 
& Ward 1973). (iii) The subsequent evolution is marked by planet's growth due to the 
accumulation of planetesimals in successive collisions. This stage is well reproduced by dy- 
namical models (Safronov 1969, Barge & Pellat 1991, 1993). When the solid core becomes 
sufficiently massive, it can accrete the surrounding gas and a giant planet, like Jupiter, is 
formed. 

However, the above scenario faces two major problems. Recent studies have shown that 
circumstellar disks are relatively turbulent and that small-scale turbulence strongly reduces 
the sedimentation of the dust particles in the ecliptic plane (Weidenschilling 1980, Cuzzi et 
al. 1993, Dubrulle et al. 1995). For particles of relevant size, the density of the dust layer is 
not sufficient to overcome the threshold imposed by Jeans instability criterion. Therefore, 

1 It can be noted that the statistical mechanics of two-dimensional turbulence is very similar to the theory 
of "violent relaxation" (Lynden-Bell 1967) by which galaxies achieve an equilibrium state. It is fascinating 
to realize that despite their very different physical nature, two-dimensional vortices and stellar systems 
share some common features. This analogy is developed in detail in Chavanis et al. (1996) and Chavanis 
(1996,1998a,1999). 
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the formation of the planetesimals, i.e the passage from cm-sized to /cm-sized particles, is 
not clearly understood. In addition, it seems difficult, with the above model, to build up 
sufficiently massive cores in less than one million years before the gas has been swept away 
by the solar wind during the T-tauri phase (Safronov 1969, Wetherill 1988). 

Both difficulties are ruled out if we allow for the presence of vortices in the disk. Their 
existence was first proposed by Von Weizacker (1944) to explain the regularity of the planet 
distribution in the solar system: the famous Titius-Bode law This idea has been reintro- 
duced recently by Barge & Sommeria (1995) and Tanga et al. (1996) who demonstrated 
that anticyclonic vortices in a rotating disk are able to capture and concentrate dust par- 
ticles. The capture is made possible by the action of the Coriolis force which pushes the 
particles inward. These results are supported by a dynamical model which integrates the 
motion of the particles in the velocity field produced by a full Navier-Stokes simulation of 
the gas component (Bracco et al. 1999, Godon & Livio 1999c). It is found that the particles 
are very efficiently captured and concentrated by the vortices. This is interesting because, 
without a confining mechanism, cm-sized bodies would rapidely fall onto the sun due to the 
inward drift associated with the velocity difference between gas and particles. Inside the 
vortices, the density of the dust cloud is increased by a large factor which is sufficient to 
trigger locally the gravitational instability and facilitate the formation of the planetesimals 
or the cores of giant planets. This trapping mechanism is quite rapid (a few rotations) and 
can reduce substantially the time scale of planet formation. 

In this paper, we present a simple analytical model for the capture of dust particles 
by coherent vortices in a rotating disk. This model is directly inspired by the numerical 
studies of Barge & Sommeria (1995) and Tanga et al. (1996) and their main results are 
recovered and confirmed. One interest of our approach is to provide analytical results 
(leading to quantitative predictions) and to isolate relevant parameters which prove to be 
particulary important in the problem. In section Q we introduce an exact solution of 
the incompressible 2D Euler equation appropriate to our sudy. This is an elliptic vortex 
with uniform vorticity matching continuously with the azimuthal Keplerian flow at large 
distances. We consider deterministic trajectories of dust particles in that vortex and derive 
analytical expressions for the capture time and the mass capture rate as a function of 
the friction parameter. In section ||, we investigate the effect of small-scale turbulence 
on the motion of the particles. Their trajectories become stochastic and their motion 
must be described in terms of diffusion equations. We estimate the diffusion coefficient 
and determine the typical length on which the particles are concentrated in the vortices. 
In section [|, we evaluate the rate of particles which diffuse away from the vortices due 
to turbulent fluctuations. An explicit expression for the "rate of escape" is obtained by 
solving a problem of quantum mechanics, namely a two-dimensional oscillator in a box. 
In appendix A, we give some details about the construction of the vortex solution and 
in appendix B, we extend Toomre instability criterion (Toomre, 1964) to the case of a 
turbulent rotating disk. 

In parallel, we apply these theoretical results to the solar nebula and make speculations 
about its actual structure. For relevant particles going from 10 cm to 100 cm in size, we 

2 An account of Von Weizsacker's theory can be found in Chandrasekhar (1946). Note that the idea of 
vortices in the solar system goes back to Descartes (1643) 
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remark that the transition between the Stokes and the Epstein regimes (at which the gas 
drag law changes) corresponds precisely to the transition between telluric (inner) and giant 
(outer) planets. Moreover, in each zone there is a prefered location where the capture 
of dust by vortices is optimal. For particles of density p s ~ 2 g/cm 3 and size ~ 30 cm (a 
typical prediction of grain growth models) , this is near the Earth orbit in the Stokes (inner) 
zone and between Jupiter and Saturn in the Epstein (outer) zone. For a broader class of 
parameters, the prefered locations cover the whole region of telluric and giant planets with a 
depletion near the asteroid belt. Inside the vortices, the surface density of the dust particles 
is increased by a factor 100 or more sufficient to trigger locally the gravitational instability. 
More precisely, we study how the surface density enhancement depends on the size of the 
particles. We show that particles must have reached at least centimetric sizes to collapse 
and form the planetesimals. Smaller particles diffuse away from the vortices on account 
of turbulent fluctuations and are not sufficiently concentrated. This implies that sticking 
processes are necessary to produce large particles. These results rehabilitate the Safronov- 
Goldreich-Ward scenario in localized regions of the disk (i.e, inside the vortices) and for 
sufficiently large (decimetric) particles . Preliminary results of this work were presented in 
Chavanis (1998b). 



2 Deterministic motion of a particle in a vortex 
2.1 The solar nebula model 

We assume that the solar nebula is disk-shaped and is in hydrostatic equilibrium in the 
vertical direction. If the mass of the disk is much less that the solar mass (< 0.1M), then 
its self-gravity can be neglected compared with the sun's attraction. In that case, the disk 
has approximately Keplerian rotation 

(D 

where r is the distance from the sun and G = 6.672 10 -8 cm 3 / (g.s 2 ) the constant of gravity. 
For thin disks, the vertical component of the solar gravity is: 

GM z . 2 
9z - ~— x ~ = ~ n z ( 2 ) 

The condition of hydrostatic equilibrium implies that the vertical pressure gradient is pre- 
cisely balanced by g z , i.e: 

() " - Pg n 2 z (3) 



dz 

is 

Pa 



If the local temperature T = j^-j^ is independant of z, then the vertical density profile of 
the gas is 



P 9 = Poe~ z2/H2 (4) 
The half-thickness (or scale height) of the disk is given by 

H = c s /Q (5) 
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where c s ~ y — is the sound speed. Other plausible temperature profiles, e.g, adiabatic 
gradient in the z direction yield similar results. 

We assume also that the gas is turbulent. Turbulence is necessary to explain the dy- 
namical evolution of the protoplanetary disk and its cooling consistent with cosmochemical 
data (Dubrulle 1993). However, its origin is not well understood and remains controver- 
sial. Various mechanisms for inducing global nebula turbulence have been proposed. Lin 
& Papaloizou (1980) and Cabot et al. ( 1987a, b) suggested that thermal convection drives 
nebula turbulence. However, the applicability of their results depends on the presence of 
abundant micron-sized dust to provide substantial thermal opacity; thus, they are ques- 
tionable when significant particle growth has already occured and the thermal opacity has 
decreased. Dubrulle (1992,1993) suggested that the Keplerian shear is the main engine of 
the turbulence. This was already pointed out by Von Weizsacker and further discussed by 
Chandrasekhar (1949) in view of the very small viscosity of the disk: "The successive rings 
of gas in the medium will have motions relative to one another, and turbulence will ensue" . 
This apparently obvious claim is actually far from straightforward to support because the 
Keplerian shear is stable with respect to linear, infinitesimal, perturbations which are usu- 
ally relied upon to induce turbulence. However, Dubrulle & Knobloch (1992) point out 
that it might be unstable to nonlinear finite amplitude perturbations, a property shared 
by most of the shear flows commonly met in engineering or laboratory experiments. The 
simplest example is the plane Couette flow, a plane parallel stream of constant shear. This 
analogy is contested by Balbus et al. (1996) who didn't evidence nonlinear instabilities in 
Keplerian disks at the respectable Reynolds numbers they achieved numerically []. These 
authors have suggested, in contrast, that turbulence in Keplerian disks could be produced 
by a powerful MHD instability (Balbus & Hawley 1991). The numerical simulations of 
Bracco et al. (1998) suggest that MHD turbulence can form magnetized vortices able 
to capture charged particles. However, in the case of the disk that is supposed to have 
spawned the solar system, it is thought that the matter was too cool to be ionized. The 
recourse to magnetic effects to render the disk turbulent is therefore suspect and the prob- 
lem of whether Keplerian disks are turbulent or not remains open. Recently, Lovelace et al. 
(1999) discoved a linear nonaxisymmetric instability in a thin Keplerian disk which can lead 
to the formation of Rossby vortices. This may open new perspectives for hydrodynamical 
turbulence in accretion disks. 

In any case, the solar nebula must have been turbulent at least during its formation 
from the collapsing protostar, because of velocity discontinuities as the infalling matter 
struck the disk. The infall probably did not stop suddenly but decayed over some interval. 
Therefore, the initial conditions in the disk were turbulent and this is sufficient to form 
vortices that survive for many rotation periods of the disk (Bracco et al. 1998,1999; Godon 
& Livio 1999a, b,c). The appearance of large-scale coherent vortices in rotating flows is 
due to the presence of the Coriolis force and the bimodal nature of turbulence (Dubrulle 
& Valdettaro 1992). At small scales, the influence of the rotation is negligible and the 
turbulence is homogeneous and isotropic. Energy cascades towards smaller and smaller 
scales up to the dissipation length. Turbulent diffusion is important and can accelerate 

3 See also the 2D simulations by Godon & Livio (1999a) starting from finite perturbations who showed 
that turbulence is not self-sustained. 
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the mixing of dust particles. At larger scales, the Coriolis force becomes dominant and the 
gas dynamics is quasi two-dimensional. On these scales, the energy transfers are inhibited 
or even reversed; this is a manifestation of the celebrated "inverse cascade" process (see, 
e.g, Kraichnan & Montgomery 1980) in two-dimensional turbulence. There is therefore the 
possibility of formation and maintenance of coherent vortices (Mc Williams 1990). These 
vortices can form either by a large-scale instability (Frisch et al. 1987, Dubrulle & Frisch 
1991, Kitchatinov et al. 1994) or by the succesive mergings of like-sign vortices, as observed 
in the simulations of Bracco et al. (1998,1999) and Godon & Livio (1999a, b,c). Due to 
the strong Keplerian shear, only anticy clonic vortices can survive this process. Initially, 
the vortices have size R <C H and typical vorticity O. Their velocity v ~ £IR is less that 
the sound speed c s = VtH and the flow can be considered as incompressible. The merging 
ends up when the Mach number M a = — reaches unity, i.e R ~ H. Bigger vortices radiate 
density waves and do not maintain (see Barge & Sommeria 1995). 

We expect therefore that, after some evolution, the disk be filled with anticyclonic 
vortices of typical size H, the disk thickness. Vortices are expected to form throughout the 
nebula and there is no reason, a priori, to beleive that certain regions of the disk should be 
excluded. However, two vortices at comparable distance from the sun will approach each 
other (due to differential rotation) and finally merge. Therefore, we do not expect more than 
one (or a few) vortices at each radial distance. On the other hand, two successive vortices 
should be separated by a distance comparable to their size R. Since R ~ H and H scales 
like a power law of the distance to the sun (r 5 / 4 in the standard model considered below), 
the distribution of vortices should be consistant with an approximate geometric progression 
of the planetary positions (Barge & Sommeria 1995). As elucidated by Graner &; Dubrulle 
(1994), the Titius-Bode law reflects more the general properties of scale invariance in the 
solar nebula than any particular physical process. 



2.2 The vortex model 

We consider the motion of a dust particle in a vortex located at a distance ro from the sun. 
For convenience, we work in a frame of reference rotating with constant angular velocity 
£1 = O(ro) and we denote by u(r,i) the velocity field of the gas in that frame. The dust 
particle is subjected to the attraction of the sun — and to a friction with the gas that 
we write — £(v — u(r, t)) where v = % is the particle velocity. We shall come back to this 



expression and to the value of the friction coefficient £ in section 2.5. Since we work in a 
rotating frame, apparent forces arise in the system. These are the Coriolis force — 2f2 A v 
and the centrifugal force Q 2 r. All things considered, the particle equation writes: 

d 2 r Jdr . ,\ nn dr / , GM\ 

This is an ordinary second order differential equation coupled to the gas motion. The case 
of a time dependant velocity field u(r, t) produced by the Navier-Stokes simulation of a 
random initial vorticity field superposed on the Keplerian rotation has been investigated 
by Bracco et al. (1999) and Godon &: Livio (1999c) who observed the capture of the dust 
particles by anticyclonic vortices. The case of a static velocity field u(r) was first considered 
by Barge & Sommeria (1995) and Tanga et al. (1996) with different vortex profiles. Barge 
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& Sommeria (1995) assume that the velocity field inside the vortex is made of concentric 
epicycles while it matches continuously the Keplerian flow at large distances. This epicyclic 
model is motivated by the underlying idea that the particles, when sufficiently concentrated, 
will retroact on the vortex and will force the gas to follow their natural motion. Tanga et al. 
(1996) consider a more complex streamfunction describing an ensemble of small vortices 
corresponding to Rossby waves corotating with the flow. This velocity field is obtained 
as a self-similar solution of the linearized equations of motion governing the large-scale 
dynamics of a turbulent nebula. These two models lead to qualitatively similar results 
indicating that dust particles are efficiently trapped into coherent anticyclonic vortices. 
However, the models differ in the details: in Tanga et al. (1996), the particles sink to the 
center of the vortices with no limit while in Barge & Sommeria (1995), the spiralling motion 
ends up on an epicycle. In that case, the friction force cancels out and the epicyclic motion 
is an exact solution of the particle equation (^). 

We must note, however, that the vortex constructed by Barge & Sommeria (1995) 
is relatively ad hoc and does not satisfy the fluid equations rigorously. In this article, 
we introduce an exact vortex solution of the incompressible 2D Euler equation and study 
analytically the motion of dust particles in that vortex. Since the vortices of the solar nebula 
are small compared with the radial distance tq (we have typically R/tq ~ 0.04, see formula 
(@D) the last term in equation (|6|) can be expanded to first order in the displacement r — tq. 
This is the so-called "epicyclic approximation" . Introducing a set of cartesian coordinates 
(x, y) centered on the vortex and such that the y-axis points in the direction opposite to 
the sun, the particle equation @ reduces to: 

At sufficiently large distances from the vortex, the velocity field is a simple shear: 

3 

u x = -My (9) 

Uy = (10) 

obtained as a first order expansion of the Keplerian velocity around ro- Its vorticity is 
u K = d x u y - d y u x = -|0. 

It remains now to specify the velocity field inside the vortex. We can construct an exact 
solution of the incompressible 2D Euler equation by taking an elliptic patch of uniform 
vorticity u. This solution is well-known by model builders (see, e.g, Saffman 1992) ^ but, 
to our knowledge, it has never been applied in an astrophysical context. Therefore, we give 
a short description of its construction in Appendix A. In the vortex, the velocity field writes 

u x = - q 2 uy (11) 

4 This solution was indicated to me by J. Sommeria. 
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where q = a/b is the aspect ratio of the elliptic patch (a, b are the semi-axis in the x and y 
directions respectively). Outside the vortex, the velocity field is given by equation ( |174|) of 
Appendix A. At large distances, we recover the Keplerian shear (p|) (|l0|) . 

The matching conditions between the elliptic vortex and the Keplerian shear require 
that u, ujk and q be related according to (see Appendix A): 

^ = 4^ (13) 
ui l + q 2 ; 

The solution (ll)fll2|) is valid for q > 1, implying < ujk/uj < 1. The vortex is anticy clonic 
(uj < 0) and is oriented with its major axis parallel to the shear streamlines. It can be shown 
to be stable to infinitesimal two-dimensional disturbances. For q — ► 1 (circular vortices), 
uj — > — oo and for q — > oo (infinitely elongated vortices), uj — ► uk- Therefore, uj is in the 
range ] — oo, — In a rotating disk, we expect that u> ~ —2Q, corresponding to a Rossby 
number of order 1. This value is achieved by vortices with aspect ratio q ~ 4, in good 
agreement with the model of Tanga et al. (1996). The epicyclic vortex considered by Barge 
& Sommeria (1995) corresponds to q = 2 and u = — Ifi. These values do not satisfy the 
matching condition (|i~3| ) with the Keplerian shear so the epicyclic vortex is not an exact 
solution of the incompressible Euler equation. It may be used, however, as an approximate 
solution if we take into account a coupling between the particles and the gas in the spirit 
of a two-fluid model. 



2.3 The capture time 

Before going into mathematical details, we recall the argument of Tanga et al. (1996) 
which shows very simply why particles are trapped by anticyclonic vortices in a rotating 
disk. Introducing a system of polar coordinates (r, 8) whose origin coincides with the vortex 
center, the radial component of equation @ reads: 

d 2 r fd0\ 2 n dO J dr \ ^ 2 2 „ 

^ = r U) +fflr S -«( S -«r)+WW« (14) 

where we have used the epicyclic approximation. The first term is a centrifugal force due 
to the rotation of the vortex (not to be confused with the centrifugal force 2 r due to the 
rotation of the disk) and the second term is the Coriolis force. The drag term in equation 
@) forces the particle velocity to approach the fluid velocity, i.e % ~ uj. For cyclonic 
vortices (uj > 0), both centrifugal and Coriolis forces are positive and push the particles 
outward. For anticyclonic vortices (uj < 0), the Coriolis force pushes the particles inward 
and comes in conflict with the centrifugal force which is always directed outward. If the 
vortex rotates rapidely, the centrifugal force prevails over the Coriolis force and the particle 
is expelled. In the other case, the Coriolis term dominates and the particle is captured. We 
conclude that only slowly rotating anticyclonic vortices can capture dust particles. This 
trapping process is specific to rotating fluids; in ordinary simulations of two-dimensional 
turbulence (without Coriolis force) the particles never penetrate the vortices. Note that the 
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epicyclic acceleration (last term in equation is always directed outward. This term 

is responsible for a flux of particles toward the sun (if the gravitational force dominates, 
i.e for y < 0) or toward the outer nebula (if the centrifugal force dominates, i.e y > 0). 
However, this flux is quite small and doesn't affect the overall trapping process (see Tanga 
et al. 1996). 

We give now a more precise analysis of the trapping process and try to derive an explicit 
expression for the capture time of the particles as a function of their friction parameter £. 
To that purpose, we notice that equations (]?]) @ with the velocity field ( |il] ) ( |l2| ) form a 
linear system of coupled differential equations. We seek therefore a solution of the form: 



Xe 



xt 



(15) 

y = Ye M (16) 

where X, Y and A are complex numbers. Substituting into equations ([?]) (||), we obtain a 
linear system of algebraic equations: 



AC A -K)A: + 2Q.( |73Y? + A ) y - «> 



3ft 



, , — -£ + 2ftA]X + (A£ + A 2 -3ft 2 )y = 
2 <?(<?-!) / 



(17) 



(18) 



There are nontrivial solutions only if the determinant of this system is zero. We are led 
therefore to a fourth order polynomial equation in A: 



A 4 + 2£A d + + ft 2 )A 2 + 3 



1 + 



q(q - 1 
1 



-£ft 2 A 



(19) 



By definition, we will say that a particle is light or heavy whether £ > ft or £ < 17 respec- 
tively. This terminology will take more sense in the sequel (see in particular section |2] 
We now consider the asymptotic limits £ — > oo and £ — » of equation (19). 
For £ — > oo (light particles), the four roots of equation (fL9j ) are 



A 



A = -£ 
317 



(double root) 

3ft 2 (g-2)(2g + l) 



(20) 



(21) 



"2(g-l) 4£ <?(g-l) 2 

The first solution is rapidely damped and will not be considered anymore. The second 
solution describes the rotation of the dust particles in the vortex: 



x = A cos 




-t/t 



capt 



t/tcapt 



(22) 
(23) 
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HQ 



In 



The particles follow ellipses of aspect ratio q and move with angular velocity — 2( q -i) 
fact, for £ — ► do, the drag term in equation (^) implies 4r — u so, in a first approximation, 
the particles just follow the vortex streamlines (their angular velocity coincides with the 
angular velocity of the fluid particles, see Appendix A). In addition, they experience a drift 
toward the center due to the combined effect of the friction and the Coriolis force. They 
reach the vortex center in a typical time: 



t 



capt 



4£ q(q ~ l) 2 
3ft 2 (q - 2)(2q + 1) 



(24) 



defined as the capture time. Note that for light particles, t cap t increases linearly with £. 
For £ — > (heavy particles), the four roots of equation (19) are: 



A 



31 + g±^/TT2^ 



A = ±tti 



q(q-l) 
( 9 -3)(2<z+l) 



(25) 



(26) 



2q(q-l) 

The second solution describes again a damped rotation of the particles but with a different 
trajectory: 

x = Acos(-Qt)e~ t/tcapt (27) 



y = — sin(-fti)e */* CBJ,< 



(28) 



The particles follow ellipses with aspect ratio 2 and move with angular velocity —ft. This is 
natural since heavy particles have the tendency to decouple from the gas and reach a pure 
epicyclic motion. However, due to a slight friction, they sink progressively in the vortex 
with a characteristic time 

t t - 2qiq - 1] - (29) 
For heavy particles, the capture time increases like £ . Very heavy particles can even 



leave the vortex. Formally, this possibility is taken into account in the first solution (25) 
which corresponds to open trajectories. The motion of heavy particles is therefore more 
complicated and demands a numerical integration of the particle equation inside and outside 
the vortex. This study will not be undertaken in that paper. We shall only consider the 
case of closed trajectories described by equations (g7|) (p8|). 

For intermediate values of £, the capture time and the angular velocity of the particles 
are plotted on figures [l] and |2] (for the particular value q = 4) . We see that t cap t presents 
an optimum at £ = £ op f. Moreover, the asymptotic expressions ( |2"4| ) and (2£) agree reason- 



ably well with the exact solution for all the values of the friction parameter. Therefore, 
considering the intersection of the asymptotes, we find that the capture time is minimum 
for 

3(9 - 2) 



Copt — 



2(g-l)(g-3) 



1/2 



ft 



and we have: 



jmm 
^capt 



3(( ? -3)(2g + l) 2 ( (? -2) 



1/2 



ft 



(30) 



(31) 
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Figure 1: Plot of t mp t vs £ for vortices with aspect ratio q = 4. The dash lines correspond 
to formulae (p3) and (29) valid for light (£ — > oo) and heavy (£ — ► 0) particles. 
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Figure 2: Plot of VL P (angular velocity of the particles) vs £ for vortices with aspect ratio 
q = 4 . Light particles (£ — > oo) move with the vortex velocity — 2 (q_i) aR d heavy particles 
(£ — ► 0) with the epicyclic velocity — $7. 
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According to equations (24) and (|29|), the condition for particle trapping (t cap t > 0) 
requires that q > 3. This implies that the vorticity must be in the range [see equation 

- -n < u) < (32) 

Particles are ejected from rapidely rotating anticyclones in agreement with the qualitative 
discussion of Tanga et al. (1996) recalled at the begining of this section. In the following, 
we shall specialize on the value of q which minimizes the capture time t™3. We find q ~ 4. 
As noted already, this value corresponds to u ~ — 2J7, i.e a Rossby number of order 1. For 
these vortices, the optimal friction parameter is t^ pt = O, and the corresponding capture 
time t™™ t = ^ is of the order of one rotation period. For £ > fi, we have 



tcapt ^ ^2 ( ligllt P articles ) ( 33 ) 



and for £ < Q 



tcapt ^ ttt (heavy particles) (34) 
3§ 

Light particles move with angular velocity — (the vortex velocity) and heavy particles 
move with angular velocity — £1 (the epicyclic velocity). These results should be compared 
with the settling time of the dust particles in the equatorial plane (see, e.g, Nakagawa et al. 
1986). Light particles are settled exponentially with a characteristic time £,/fl 2 comparable 
with equation (33). On the other hand, heavy particles undergo overdamped oscillations 
around the midplane with a period of the order of O -1 and a settling time 2/£ similar to 
expression (34). We shall come back on the analogy between the vortex trapping and the 



dust settling in section 3.1 



The models of Barge & Sommeria (1995) and Tanga et al. (1996) lead to qualitatively 
similar results. In the model of Tanga et al. (1995), the particles sink to the center of 
the vortex on a typical time t cap t{t]) which also presents a minimum when £ is of order 
n. Moreover, for light particles, t cap t increases linearly with £ like in equation (|33|) . The 
model of Barge & Sommeria (1995) is physically different since the particles do not really 
reach the vortex center but end up on an epicycle after a time ~ l/£. However, the general 
tendancy is the same. Light particles (£ 3> fi) mainly follow the streamlines of the gas 
and stop on epicycles close to the vortex edge (like in case (a) of their figure 1). Optimal 
particles with friction parameter £ ~ Q reach deeper epicycles, almost at the center of the 
vortex (case (b)). Heavy particles (£ <C f2) take a long time to connect an epicycle and can 
even escape from the vortex since their motion is nearly unaffected by the friction drag. 
This is also a possibility in our model and in Tanga et al. (1996). Therefore, the three 
vortex models give relatively similar results even if their physical contents are different. 
The recent numerical simulations of Godon & Livio (1999c) for a compressible flow also 
show a capture optimum when the drag parameter is of the order of the orbital frequency. 

2.4 The mass capture rate 

In addition to the capture time, an important aspect of the problem concerns the mass 
capture rate (Barge & Sommeria 1995). This quantity can be estimated as follows. First, 
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we have to determine the capture cross section of the vortex as a function of the friction 
parameter £. Without the effect of the drift, a particle with impact parameter r\ would cross 
the vortex in a typical time t cross ~ R/uk, where uk = 5^77 is the Keplerian velocity. The 
particle will be captured by the vortex if, during that time, the deflection due to the drift 
is precisely of order rj. Light particles (£ — > 00) follow the edge of the vortex and, for them, 
Vdrift ~ R/tcapt- On the other hand, heavy particles (£ — > 0) keep their rectilinear motion 
and enter directly the vortex so that, for them, v drift ~ v/tcapt- The critical parameter r/ c 
is given by the condition t cross x v drift ~ Moreover, for optimal particles with £ ~ fi, 
we expect that r/ c ~ i?. Regrouping all these results, we obtain 



1/2 

/(£) = g - ( j J (light particles) (35) 

/(0 = ^ - ^ (heavy particles) (36) 

These results agree with the asymptotic behaviour of the function /(£) determined numer- 
ically by Barge & Sommeria (1995) in their model. 

The mass capture rate can be estimated by assuming that the particles are carried to 
the vortex by the Keplerian shear uk = and that they are continuously renewed 

by the inward drift (directed towards the sun) associated with the velocity difference AV 
between gas and particles (Barge k, Sommeria 1995). Another alternative, consistent with 
the simulations of Bracco et al. (1999), is that the particles are already concentrated 
by the turbulent fluctuations that accompany vortex formation in the early stages of the 
protoplanetary nebula. It is likely that both mechanisms come into play in the capture 
process. Considering the first possibility, which can be studied analytically, we have (Barge 
& Sommeria 1995): 

^ = 2 ^ a d u K dy = ±<T d nR 2 f 2 (0 (37) 

where ad is the surface density of the dust particles outside the vortex. The total mass 
collected by the vortex during its lifetime is 

M life = ^(nt life )a d R 2 f 2 (0 (38) 

The mass capture rate is proportional to the effective surface ~ / 2 (£)i? 2 of the vortex. As 
first emphasized by Barge &i Sommeria (1995), it is maximum for particles with £ ~ Q. 

In conclusion, various vortex models and different physical arguments show that the 
capture is optimal for particles whose friction parameter is close to the disk rotation. We 
shall now consider the implications of this result on the structure of the solar nebula. 

2.5 Application to the solar nebula 

We shall assume that the solid particles are spherical, of radius a and density p s . The value 
of the friction parameter £ depends whether the size of the particles is larger or smaller 
than the mean free path A in the gas (see, e.g, Weidenschilling 1977). When a < |A, we 
are in the Epstein regime and: 

* = P- (39) 
2p s a 
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where a g is the gas surface density. 

When a > |A, the situation is more complicated because, in general, the friction pa- 
rameter depends on the velocity difference |v — u| between the particles and the gas. There 
is, however, a regime where ^ is independant of the particle relative velocity. This is the 
Stokes regime in which: 

This regime is valid when the particle Reynolds number R p = jjr|v — u|, where v g = ^c s \ 
is the kinematic viscosity of the gas, is smaller than 1. This is the case for light particles 
which move typically with the gas velocity v ~ u. For heavy particles on the contrary, 
we have R p > 1 and, rigorously speaking, the friction parameter depends on the relative 
velocity. The use of a more exact expression for £ would require a numerical integration 
of the particle trajectory but the results shouldn't dramatically differ from those obtained 
with expression (^0|). Since we are mainly interested by orders of magnitude, we will use 



expression ( |40| ) for all particle sizes (with a > |A), but keep in mind this uncertainty for 
large particles. 

We shall adopt a standard model of the solar nebula, following Cuzzi et al. (1993). It 
corresponds to a "minimum mass" circumstellar nebula with parameters: 

/ \ -3/2 

n = 27T \TMj) years_1 (41) 

ad = l \ihj) g/cm2 (42) 

<T a = Vm(— -) g/cm 2 (43) 



1A.U 



H = 0.04 



r \ 



5/4 



IA 



— ) A.U (44) 



11/4 



x = \Tnj) cm <45) 

We take 1 year = 3.15 10 7 s and 1 A.U = 1.49 10 13 cm. 

For a given type of particles, the transition between the Epstein and the Stokes regime 
is achieved at a specific distance from the sun given by: 

j_ 

V'a.U (46) 
9 1cm y v ' 

We are particularly interested by particles of order 10 cm in size. Indeed, smaller particles 
can grow by aggregation processes without the aid of vortices. However, when they reach 
decimetric sizes, collisional fragmentation becomes prohibitive and prevents further evolu- 
tion (on relevant time scales) . It is therefore at this range of sizes that the vortex scenario 
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should come at work and facilitate the formation of bigger structures, the so-called plan- 
etesimals. For particles between 10 and 100 cm in size, we find that the critical radius ( |46| ) 
at which the gas drag law changes lies in the range: 

1.7A.U <r c < 3.9A.U (47) 

that is to say just at the separation between telluric (inner) and giant (outer) planets. This 
result is relatively robust because it depends only on a small power of the particles size 
and is independant on their density. We can therefore wonder if there is not a connection 
between the division of the solar system in two groups of planets (telluric and gaseous) and 
the two regimes (Stokes and Epstein) of the gas drag law in the primordial nebula. In the 
following we show how the vortex scenario can give further support to this idea. 

Since the friction coefficient of a particle with size a and density p s depends on param- 
eters (like o~ g , A and f2) which are functions of the distance r from the sun, the friction 
coefficient £ is itself a function of r. Combining equations ((35^) (|40|) with equations (fH] ) 
OH) OH), we obtain: 

J- = — ~ — r* if r < r c (Stokes zone) (48) 
U a z p s 

— = r 2 if r > r c (Epstein zone) (49) 

il ap s 

where r is measured in A.U, a in cm and p s in g/cm 3 . 



According to the results of section |2.3| , the capture time ilt cap t (measured in rotation 
periods) is optimal when £/il = 1. Since is a function of r with a maximum at r c , this 
criterion determines two prefered locations in the disk, one in each zone (see figure |3|). In 
the Stokes (inner) zone, the optimum is at: 



and in the Epstein (outer) zone, it is at: 



2 \ 4/5 

a Ps N 
1913 



/850\ 2/3 



(50) 



r ut =[ (51) 



More generally, we can combine equations (p3|) (p3|) with equations (fig) (f49j) to express the 
capture time ilt cap t as a function of the distance to the sun (for a given type of particles). 
We find a W-shaped curve (see figure |j) determined by the power laws: 

Sitcom = ^r" 5 / 4 (r < r m ) (52) 

ttt capt = ^i r 5 / 4 (r in <r <r c ) (53) 
a 2 p s 

ttt capt = ' 2 ^L r - z l 2 (r c <r < r out ) (54) 
ap s 
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Figure 3: Evolution of the friction parameter £/f2 as a function of the distance to the sun 
for a given type of particles (the figure corresponds to particles with size a = 30 cm and 
bulk density p s = 2 g/cm 3 ). The friction parameter is maximum at r = r c where the gas 
drag law passes from the Stokes to the Epstein regime. The condition = 1 determines 
two optimal regions in the disk. 
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Figure 4: Evolution of the capture time £lt cap t as a function of the distance to the sun for 
a given type of particles (a = 30 cm, p s = 2 g/cm 3 ). We have represented the present 
position of the planets. The capture time is optimum near the Earth (in the Stokes zone) 
and between Jupiter and Saturn (in the Epstein zone). 
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Table 1: Minimum sizes of dust particles (in cm) which can trigger the gravitational insta- 
bility (p s = 2g/cm 3 ). 



Flanets 


1 A TT\ 

r(A.U ) 


A 




„ cone 


capt 
n ■ 
mm 


O'min 


Telluric planets : 














Mercury 


0.387 


0.07 


17 


14 


7 


3 


Venus 


0.723 


0.4 


25 


20 


11 


5 


Earth 


1 


1 


31 


25 


13 


6 


Mars 


1.52 


3 


40 


33 


17 


8 


Asteroid belt : 


3 


20 
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50 


15 
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Giant planets : 














Jupiter 


5.20 


93 


36 


23 


7 


1 


Saturn 


9.52 


492 


14 


9 


3 


0.5 


Uranus 


19.2 


3364 


5 


3 


1 


0.2 


Neptune 


30.0 


11523 


3 


2 


0.5 


0.1 



nt capt = |gr 3 / 2 (r > r out ) (55) 

In order to make a numerical application, we assume that all the particles have the 
same density p s ~ 2g/cm 3 (the density of a composite rock- ice material) and the same size 
a ~ 30cm (a typical prediction of grain growth models). In principle, we should consider a 
spectrum of sizes and densities but we choose these particular values in order to illustrate 
at best the predictions of the vortex model. For these values, the optimum in the inner 
zone is at rj n ~ 1A.U, i.e near the Earth orbit, and the optimum in the outer zone is at 
T out ~ 6A.U, i.e between Jupiter and Saturn's orbits. The transition between the Stokes 
and the Epstein regime occurs at r c ~ 2.7A.U, i.e near the asteroid belt. 

Alternatively, we can determine, at each heliocentric distance, the size of the particles 
which are preferentially concentrated by the vortices. The results are indicated on table [l] 
(see also figures || and They show that the optimal sizes lie between 1 and 50 cm in 
the region of the planets (for a bulk density p s = 2g/cm 3 ). Such particles are concentrated 
in the vortices after only one rotation period. By contrast, the capture time for particles 
of lpm in size (the initial size of the dust grains in the primordial nebula) is Qt cap t ~ 10 9 
(we have a similar characteristic value for their settling time). This exceeds the lifetime of 



a circumstellar disk by many orders of magnitude. These results (see also sections 3.4 and 



Q|) indicate that sticking of particles up to cm-sized bodies is an indispensable step in the 
process of planet formation. 

We can also study how the mass captured by the vortices varies throughout the nebula 
and how it depends on the size of the particles. According to equations ( |38[ ) ([42]) and (|44]), 
the mass captured by a vortex during its lifetime can be written: 

= 8.916 W- 4 (nt Ufe )rf(0 (56) 
where we have introduced the Earth mass M ffi = 5.976 10 27 g as a normalization factor. 
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Figure 5: Variation of the friction parameter with the size of the particles at r = 1 A.U 
and r = 5 A.U (we have taken p s = 2g/cm 3 ). The capture is optimum for particles with 
friction parameter £ ~ ft. This corresponds to decimetric sizes in the region of the planets. 
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Figure 6: Variation of the capture time with the size of the the particles at r = 1 A.U 
and r = 5 A.U (p s = 2g/cm 3 ). The capture time is minimum for particles with radius 
a ~ 30cm. 
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Figure 8: Variation of Mnf e with the heliocentric distance for particles with size a = 30cm 
and bulk density p s = 2g/cm 3 . 
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At r = 1A.U and for optimal particles with size a ~ 30cm and bulk density p s ~ 2g/cm? 
(see table |l]), we have £ ~ £1 leading to rf 2 ~ 1. Since the vortex lifetime is not accurately 
known, it makes sense to calibrate its value so as to satisfy Muf e (lA.U) = M®. This yields 
Qtufe — 1122, corresponding to ~ 180 rotation periods, in agreement with the estimate of 
Barge & Sommeria (1995) based on the value of the disk a-viscosity and with the numerical 
simulations of Godon & Livio (1999a,b,c). Figure ^ shows how the mass Muf e depends on 
the size of the particles at a given position of the solar nebula. We can also determine how 
it varies with the heliocentric distance for a given type of particles. The results are reported 
on figure || (full line) for particles with a ~ 30cm and p s ~ 2g/cm 3 . The captured mass 
presents a global maximum near Jupiter's orbit and a plateau in the Earth's region. These 
results complete the work of Barge & Sommeria (1995) who first noticed the existence 
of an optimum near Jupiter. However, they didn't consider the Stokes regime in their 
article and made the wrong statement that "inside Jupiter's orbit, particle concentration 
occurs in an annular region at the vortex periphery [because the friction parameter £ is 
larger]". In reality, the friction parameter decreases anew when we enter the Stokes domain 
at r < r c (see figure ||). This allows the existence of another optimum near the Earth 
orbit. Therefore, the mass collected by the vortices in the inner zone is larger that one 
would obtain without this change of regime (the dash line in figure ||| would be found if the 
Epstein law was used in the whole disk). This transition may explain the division of the 
solar system in two groups of planets. Moreover, the vortex scenario explains naturally the 
disymmetry between these two groups: the mass capture rate is larger in the outer zone 
simply because the vortices are larger. Indeed, the captured mass is not only proportional 
to f 2 (which would give the symmetrical dotted line of figure ||) but also to the product 
(JdR 2 which increases linearly with r. This effect may explain the difference in size and 
mass between telluric and giant planets. Moreover, the mass captured by the vortices in 
the outer zone should be larger than these estimates since they intercept in priority the 
matter drifting towards the sun. On the other hand, the intermediate region at ~ 3 A.U 
should be further depleted due to its proximity with the global maximum. 

In conclusion, we find that there are two locations in the primordial nebula where 
the trapping of dust by vortices is optimal. These locations belong to the Stokes and to 
the Epstein zones are fall near the Earth and Jupiter's orbit respectively. The zone of 
transition of the gas drag law is consistent with the position of the asteroid belt which 
marks the separation between telluric and giant planets. The asymmetry between the two 
groups of planets may be related to the size of the vortices which are bigger in the outer 
zone and therefore capture more mass. The exact values of rj n , r c and r ou t depend on the 
spectrum of size of the particles, which is not well-known, but the W-shape of figures |I| and 
|H is generic and agrees with the global structure of the solar system. 

3 Stochastic motion of a particle in a vortex 
3.1 The diffusion equation 

Due to small-scale turbulence, the motion of a particle in a vortex is not deterministic 
but stochastic. Turbulent fluctuations produce some kicks which progressively deviate the 
particle from its unperturbed trajectory. This is similar to what happens to a colloidal 
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particle in suspension in a liquid (Brownian motion) . An individual fluctuation has a minute 
effect on the motion of the particle, but the repeated action of these fluctuations gives rise 
to a macroscopic process of diffusion. The effect of turbulence on the sedimentation of 
dust particles in the protoplanetary nebula has been considered in detail by Dubrulle et al. 
(1995). We use a similar approach to study the effect of turbulence on the capture of dust 
by large-scale vortices. 

The transport equation governing the evolution of the dust surface density inside a 
vortex can be written: 

^ + V(a d (v))=V(DVa d ) (57) 
where (v) is the mean velocity of the particles and D their diffusivity. The mean velocity 



(v) is given by the deterministic model of section 2J3. According to equations (|22|) ([23|) or 
(|^) (|D, it can be written: 

(v)=V-J-r (58) 

leapt 

where V corresponds to the pure rotation of the particles in the vortex and —r/t cap t is their 
drift towards the center. In the case of light particles (£ 3> 0), V is equal to the velocity of 
the vortex while in the case of heavy particles (£ <C f2), V is equal to the epicyclic velocity 



(see section 2.3). When (|5q ) is substituted into (|57|), the diffusion equation takes the form: 



^ + V(a d V) = V ( DVa d + ^-v) (59) 

Ot \ t ca pt J 



The first term in the r.h.s is a pure diffusion due to small-scale turbulence and the second 
term is a drift toward the vortex center due to the combined effect of the Coriolis force 
and the vortex (anticyclonic) rotation. Equation ( |59] ) illustrates the bimodal nature of 
turbulence: the diffusion is due to small-scale fluctuations hardly affected by the rotation 
of the disk (3D turbulence) and the trapping process is a consequence of the Coriolis force 
and the existence of coherent structures in the disk (2D turbulence). 

Since we are mainly interested by orders of magnitude and in order to avoid unnecessary 
mathematical complications, we shall consider from now on that the vortices are circular 
with typical radius R. In this approximation, we can restrict ourselves to axisymmetric 
solutions for which the advection term in equation (^) cancels out. We are led therefore 
to study the Fokker-Planck equation: 

d ° d v( D Va d +^Lr) (60) 

capt J 



dt \ t 



For an initial condition consisting of a Dirac function centered at ro, there is a well-known 
analytical solution of equation (^) : 

where M is the total mass of particles contained in the vortex. This formula shows that 
the relaxation time is equal to the capture time t cap t not affected by turbulence. 
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The equilibrium solution of the Fokker-Planck equation (|60|) satisfies the condition: 

DVa d + -^-r = (62) 



t 



capt 



expressing the balance between the diffusion and the drift. It corresponds to a gaussian 
density profile of the form: 



„2 



a d = a d (0)e 2a W (63) 

This is also the solution of equation ( |6ll) for t —> +oo. In practice, the equilibrium distri- 
bution ( |63[ ) is established for t ~ t C apt- Then, the particles are concentrated in the vortices 
on a typical length: 

Id = V 2Dt capt (64) 

3.2 Vertical sedimentation 

Equation (|S0D is similar to the diffusion equation 

dn d d f dn d fl 2 \ 



dt dz\dz £ 

used by Dubrulle et al. (1995) to describe the sedimentation of the dust particles and 
determine the sub-disk scale height (see also Weidenschilling 1980). The drift towards the 
vortex center is replaced in their study by the drift — ^-z towards the ecliptic plane due to 

gravity. For light particles, t cap t ~ Jf? [see equation (^3|)] and the two equations coincide 
up to numerical factors. This implies that the particles are concentrated in the vortices 
on a lenght l d comparable with the sub-disk scale height H d determined by Dubrulle et al. 
(1995) [see their section 3]. 

It is relatively straightforward to include the vertical sedimentation of particles in our 
study, although we shall specialize in the following on their horizontal accumulation in 
vortices. Introducing the volume density p d (r, z, t) = -^n^z, t)o~d(r, t) and using equations 
(|60|) and (pq) , we obtain: 

01,11 v(l»V P , + ^p d ( lr x + *)\ (66) 



dt I £ 

where and z are the component of r in the directions perpendicular and parallel to the 
disk rotation vector Q. Integrating equation ( p6| ) on the vertical direction returns equation 
(|60|) for the surface density. 



3.3 Diffusivity of dust particles 

It remains now to specify the value of the diffusion coefficient appearing in equation (|60|). 
In general, the turbulent viscosity of the gas is written under the form v ~ a-^ where c s 
is the sound speed and a a non dimensional parameter which measures the efficiency of 
turbulence (Shakura & Sunyaev 1973). Following current nebula models, 10 -4 < a < 10 -2 . 
When the turbulence is generated by differential rotation, a = 2 10 -3 (Dubrulle 1992) and 
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we shall take this value for numerical applications. Since the disk height is H ~ c s /Q (see 
equation (|5|)), the turbulent viscosity can be written v ~ aH 2 Q or, alternatively, v ~ aflR 2 
where R is the vortex radius. More precisely, assuming a power law spectrum E(k) ~ A: -7 
for the gas turbulence (with 7 ~ 5/3) , we have 

°°* (67) 



According to Dubrulle et al. (1995), the diffusivity of the particles can be written: 

D = g{g)v (68) 

where 

+ B k0 J { ' 

is a function of the friction parameter. The reduction factor of Volk et al. (1980): 

- £fc (70) 

depends on the size k^ 1 ~ \fa¥l of the largest eddies of turbulence and on the systematic 
velocity v s of the dust grains. In the vortices, v s is equal to the drift velocity r/t capt . 
For light particles (£ 3> fi), g(£) — > 1 and 

D = (light particles) (71) 

V7+ 1 

This result is expectable since light particles mainly follow the streamlines of the gas. 
Consequently, their diffusivity tends to the gas turbulent viscosity v. 
For heavy particles (£ <C £1), g(£) — > y£A^ an d 

D = ^yrpy J (heavy particles) (72) 



For £ — > 0, I? — > since there is no coupling with the gas. For £ ~ O, equations ( |7l| ) and 
( [72] ) give the same result, so we can use these expressions in the whole range of friction 
parameters. Note, however, that the diffusion approximation is not strictly valid for heavy 
particles, so equation (72) must be taken with care. 

Substituting equations (33) (34) and (|7l|) (|72|) into equation ([64]), the concentration 
length can be written explicitely: 



ld_( 16a y (t 

r \n 



l/2/t \ V2 

(light particles) (73) 



/ i6« y/Vny/* 

y 3 ^TT J \JJ ( heav y particles) (74) 



i d ( i6a v f yn\ 

R 
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As expected, the concentration is optimum for £ ~ 0. In that case, the particles are 
distributed over a length ~ yfaR. Lighter and heavier particles are less concentrated. 
When 

* > 3 -^H ^ 153 (75) 



n 16a 

or 



2 



^I^ttJ^ 410 <76) 

the drift is negligible and there is no concentration (l^ > R). 
3.4 Application to the solar nebula 

We now apply these results to the case of the solar nebula and show how the vortex scenario 
can make possible the formation of planetesimals at certain prefered locations of the disk. 

It is generally beleived that planetesimals formed by the gravitational instability of the 
particle sublayer (Safronov 1969, Goldreich & Ward 1973). The dispersion relation for an 
infinite uniformly rotating sheet of gas is (see, e.g, Binney & Tremaine 1987): 

to 2 = k 2 c 2 d + 4Q 2 - 27rGa d k (77) 

where a d is the unperturbed surface density of the particles and c d their velocity dispersion. 
The particle sub-layer (which behaves as a very compressible fluid) will be unstable provided 
that uj 2 < 0, for some k. From equation (f77|), we obtain the criterion for gravitational 
instability (Toomre, 1964): 

Cd < = Vcrit (78) 

Once instability is triggered, the system crumbles into numerous planetesimals of order 
10 km in size. Moreover, the growth time of density perturbation is predicted to be short, 
of the order of an orbital period. In addition, the instability criterion gives the impression 
that its operation does not require any sticking mechanism. Goldreich & Ward (1973) 
state that "...the fate of planetary accretion no longer appears to hinge on the stickiness 
of the surface of dust particles". This is very attractive because sticking mechanisms are 
relatively ad hoc and ill- understood. For these reasons (and also for a lack of alternatives), 
the Safronov-Goldreich-Ward scenario was nearly universaly accepted as the key mechanism 
for forming planetesimals. Therefore, a swarm of bodies of a few km in diameter was a 
common starting point for numerical simulations of planetary formation. 

However, Weidenschilling (1980), followed by Cuzzi et al. (1993) and Dubrulle et al. 
(1995), realized that this simplisitic picture was ruled out if the primordial nebula was tur- 
bulent. Indeed, turbulence reduces considerably the vertical sedimentation of the dust 
particles and prevents gravitational instability. According to equations (|4l| ) and (]4^), 
the velocity threshold imposed by the instability criterion ( |78| ) is of order V cr u = 5cm /s 
throughout the nebula. Even if the nebula as a whole was perfectly laminar, the formation 
of a dense layer of particles (considered as a heavy fluid) would create a turbulent shear 
with the overlying gas (see, e.g., Weidenschilling & Cuzzi 1993). The velocity dispersion of 
the particles can be estimated by Cd ~ \/ac s and easily reaches several meters per second 
(when a numerical value is needed, we shall take Cd = 5m/ s). Therefore, the instability 
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criterion ( |7^ ) is not satisfied (see appendix B for a rigorous derivation of this criterion in 
a turbulent disk). Turbulence is responsible for too high velocity dispersions or, alterna- 
tively, the surface density of the dust sublayer is not sufficient to trigger the gravitational 
instability. The density needs to be increased by a factor A c ~ 100 or more to overcome 
the threshold imposed by Toomre instability criterion. 

There is therefore a major problem to form planetesimals by gravitational instability in 
a turbulent disk. As first suggested by Barge & Sommeria (1995), the presence of vortices 
in the disk can solve this problem. Indeed, by capturing and concentrating the particles, 
the vortices can increase locally the surface density of the dust sublayer and initiate the 
gravitational instability. Let us first discuss the concentration effect. Inside a vortex, an 
initial mass of order ad^R 2 is concentrated on a typical length Id given by equation pif) . 
The surface density is therefore amplified by a factor (R/ld) 2 depending on the size of the 
particles. Using equations (73)(fT4]), this amplification can be written explicitely 



a vort \ 3 

= rrz — \/7 + 1— (light particles) (79) 

J cone. 16 « f 



o-d 

/ a vort\ 3 /£\l/2 

(^LcT^^W (»«vy parti*,) (80) 
The amplification is maximum for £ ~ ft and takes the value A^x = tJ-\/7 + 1 — 150. As 



we have seen in section [2.5| , this corresponds to particles of size a = 30cm and bulk density 
Pd = 2g/cm? in the regions of the Earth and Jupiter. This enhancement is sufficient to 
satisfy the instability criterion (|78|) 0. For microscopic particles, on the contrary, there is 
no density enhancement. In that case, l d ~ 10 3 R (in the region of the planets) and the 
particles rapidely diffuse away from the vortices (see section |4~5| ) . Gravitational instability 
will be possible provided that: 

I < -L^i * 153 (81) 

On table [l], we report, as a function of the heliocentric distance, the minimum size of the 
particles which satisfy this criterion (see also figures ^ and p|). These results indicate that 
particles must have grown up to some centimeters to trigger the gravitational instability. 
Therefore, sticking processes are needed to reach this range of sizes (recall that this was 
claimed to be not necessary in the initial Safronov-Goldreich-Ward scenario). 

In conclusion, by allowing a local enhancement of the particle surface density, the vor- 
tices can favour the formation of planetesimals by gravitational instability. This rehabili- 
tates the Safronov-Goldreich-Ward theory at certain prefered locations of the disk (i.e inside 
the vortices) and for sufficiently large (decimetric) particles. A sufficient enhancement is 
achieved simply by the horizontal concentration of the dust layer in the vortex (a process 
similar to the vertical sedimentation). However, this mechanism alone is not sufficient to 
produce enough planetesimals to form the planets. The vortices must also capture the 

5 We can argue that when Id is comparable to the sub-disk thickness Hd, the instability criterion derived 
in the case of a thin disk is not applicable anymore. However, according to Jeans criterion, a volume element 
of size Hp and mass M p is unstable if the velocity dispersion of the particles c? d is less than their potential 
energy G ^ d . Since Cd ~ H d U and Md ~ o- d H\, this condition returns the criterion (F§|). 
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Figure 9: Enhancement of the dust surface density inside the vortices (due to concentration) 
as a function of the size of the particles at r = 1A.U and r = 5A.U (p s = 2g/cm 3 ). 
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Figure 10: Enhancement of the dust surface density inside the vortices (due to capture) as 
a function of the size of the particles at r = 1A.U and r = 5A.U (p s = 2g/cm 3 ). 
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surrounding mass. This is another source for density enhancement (Barge & Sommeria 
1995). Returning to equation (38), we find that the average surface density of the particles 
collected by a vortex during its lifetime is 



vort\ 



capt. 



(82) 



with f2i 



life 



1122 (see section |2.5| ). The maximum amplification, reached by particles 
with £ ~ f2, is Amax = ^{^Mife) — 536 a little bit larger than the previous value (recall 
however that tnf e is not known precisley). Gravitational instability will be possible for 
particles whose friction parameter satisfies 



£ 3 1 



5.4 



(83) 



See also table |l| and figures || |To| for the same criterion expressed in terms of the size of the 
particles. 

If we now take into account both the concentration effect and the capture process, we 
obtain an amplification 
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~ 28.6 



(85) 



but, even in this optimistic situation, the particles must have reached relatively large sizes 
to trigger the gravitational instability (see table |l]). Of course, if the vortex lifetime is 
increased, smaller particles have the possibility to collapse since the vortex captures more 
mass. In fact, this is not completely correct because the previous results assume that the 
escape of particles due to turbulent fluctuations can be neglected. This is not always the 
case (in particular for small particles) and this problem is now considered in detail. 



4 The rate of escape 

4.1 Formulation of the problem 



The diffusion equation ( |60D is similar, in structure, with the Kramers-Chandrasekhar equa- 
tion: 

9/ 9 („B} 



St o,\ D t + u V (86) 

introduced in the case of colloidal suspensions and in stellar dynamics (Kramers 1940, 
Chandrasekhar 1943a, b). In this equation, f(y,t) governs the velocity distribution of the 
particles in the system. The first term in the r.h.s is a pure diffusion and the second term 
is a dynamical friction. These terms model the encounters between stars or the collisions 
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between the colloidal particles and the fluid molecules. Comparing equations ( |60| ) and 
we see that the position in ( pCf ) plays the role of the velocity in (|86|) and the capture time 
tcapt the role of the friction time In particular the friction force and the drift term are 
linear in v and r respectively. Equation ( |S6| ) was used by Chandrasekhar (1943b) to study 
the evaporation of stars in globular clusters (this is similar to the Kramers problem for the 
escape of colloidal suspensions over a potential barrier). Due to collisions, some stars may 
acquire very high energies and escape from the system (being ultimately captured by the 
gravity of nearby objects). Similarly, in our situation, turbulent fluctuations allow some 
dust particles to diffuse towards higher and higher radii and finally leave the vortex (being 
eventually transported by the local Keplerian shear). In each case, the friction force or the 
drift acts against the diffusion and can reduce significantly the escape process. 

We try now to evaluate the rate of particles that leave the vortex on account of turbulent 
fluctuations. To that purpose, we formulate the problem in terms of the density probability 
W(r,t) = W(\r\,t) that a particle located initially in the annulus between |r| = r$ and 
tq + drQ will be found in the surface element around r at time t. According to equation 
(|60|), the time evolution of the probability W(r,t) is given by 

dW ( W \ 

V I DVW + — — r ) (87) 



dt V t. 



with initial condition 



capt 



W(v,t) = 6{lr }- ro) as t^O 

27T7-0 



where 5 stands for Dirac's 5-function. We assume that when the particle reaches the vortex 
boundary at |r| = R, it is immediately transported away be the Keplerian shear. In other 
words, we adopt the boundary condition: 

W(r,t) = for \r\ = R for all i >0 (89) 

We call 



( 



W 

DVW + r (90) 



V tcapt 

the current of probability, i.e Jdlh gives the probability that a particle crosses an element 
of length dl between t and t + dt (n is a unit vector normal to the element of length under 
consideration). 

We first introduce the probability p(ro, t)dt that a particle located initially in the annulus 
between |r| = ro and ro + dro reaches for the first time the vortex boundary between t and 
t + dt. According to what was just said concerning the interpretation of (|90|), we have: 



p(r ,t)=i Jhdl = -( ) (91) 
J\v\=R V dr J r=R 

The total probability Q(ro,t) that the particle has reached the vortex boundary between 
and t is ^ 

Q(r ,t) = f p(r ,t')dt' (92) 



o 
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Finally, we average Q(ro,t) over an appropriate range of initial positions in order to get 
the expectation Q(t) that the particle has left the vortex at time t. We have 



Q{t) = j Q(r ,t)n(r )2irr dr (93) 

where /i(|ro|) governs the initial probability distribution of the particles in the vortex. In 
terms of the function Q, the rate of escape of the particles can be written 

_L = ^_dQ_ (94) 
t esc (1 - Q) dt K ' 

As mentioned already, this problem is similar to the diffusion of colloidal suspensions 
over a potential barrier or to the evaporation of stars in globular clusters. As will soon be- 
come apparent, it reduces to solving a pseudo Schrodinger equation for a quantum oscillator 
in a "box" . In this analogy, the rate of escape appears to be related to the fundamental 
eigenvalue of the quantum problem. An explicit expression for the rate of escape can be 
obtained in two limits: when £ — » or £ — » do, the drift term can be ignored and the 
Fokker-Planck equation reduces to a pure diffusion equation (section ^^). On the other 
hand, when £ ~ 0, the drift term is dominant and a perturbation approach inspired by 
the work of Sommerfeld and Chandrasekhar can be implemented to determine the ground 
state of the artificially limited quantum oscillator (sections and f4.4| ). 

4.2 The rate of escape when the drift term is ignored 

When the drift term can be ignored (i.e for very light or very heavy particles, see inequalities 
(|7q) ([rq) ), the Fokker-Planck equation (p7|) reduces to a pure diffusion equation 

dW 

= DAW (95) 
dt K ' 



which has to be solved in a circular domain with boundary conditions ( p8|) and (|89|). The 
solution of this classical problem is 

1 + °° 1 f Tq 

nR J l (°0n) V R 

xJo{a 0n -je b2 (96) 

where J m is Bessel function of order m and the aon's denote the roots of Bessel function 
Jo- The probability p(ro, t) that a particle with an initial position r$ has reached the vortex 
boundary between t and t + dt is 

i . +00 , \ 2 

2D ^ « n T I r \ _ Da 0n t 
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and the total probability that the particle has escaped during the interval (0, i) is 

+ °° 1 / r \ 

Q(r ,t) = 2 V — -J ( a 0n -^ ) 



D,>7, 



x^l-e __ ^~J (98) 

Averaging the foregoing expression over all ro's in the range [0, R] with equal probability 
/i(ro) = 1/ttR 2 (this corresponds to an initially homogeneous distribution of particles in 
the vortex), we obtain 

W=4£^(l-e-^H (99) 

To sufficient accuracy, we can keep only the first term in the series. The expectation that 
the particle has left the vortex at time t is therefore: 



(100) 



Since aoi — 2.40482..., this term represents ~ 70% of the value of the series (^) and the 
approximation ( |100| ) is reasonable. 

In conclusion, when the drift term is ignored, we find that the escape time is 

R , 

(101) 



Da 01 



It corresponds, typically, to the time needed by the particles to diffuse over a distance ~ R, 
the vortex size. For light particles, using (fTi"|), we obtain explicitely 



VLt esc = (light particles) (102) 



and for heavy particles 



1/2 

Qt esc = v '„ — ( — ) (heavy particles) (103) 



\/tTT/o 



We shall come back to these expressions in section 4.5 



4.3 The effective Schrodinger equation 

We now return to the general problem for the rate of escape when proper allowance is made 
for the drift. We find it convenient to introduce the notations 

t = — - — and p= — f (104) 
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or, in words, the time is measured in terms of the capture time and the distances are 
normalized by the concentration length (Q). We let also: 

w = 2Dt capt W (105) 

and 

' R (106) 



y/2Dt cap t 

In terms of these new variables, the problem (|S7|)(|88D and ( p9| ) takes the form: 

dw I 

— = -Aw + V(wp) (107) 

OT 2 

w fr T ) = S j£p£*L as r ^ (108) 

w(p, t) = for p = poo for all r > (109) 
With the change of variables 

w = ^ e - p2/2 (110) 



we can transform the Fokker-Planck equation ( 107 ) into a Schrodinger equation (with 
imaginary time) for a quantum oscillator: 

^ = ^ + (1-1^ (111) 

However, contrary to the standard quantum problem, equatio n (|111| ) has to be solved in 
a bounded domain of size p^ with the boundary condition ( |109|) , In other words, our 
problem consists in determining the characteristic functions of a quantum oscillator in a 
"box" . 

First, we notice that a separation of the variables can be effected by the substitution 

tf> = <P(p)e- Xr (112) 

where A is, for the moment, an unspecified constant. This transformation reduces the 
Schrodinger equation in a second order ordinary differential equation: 

d 2 6 1 d6 , , 9 . , 

J + - P i +(2+2X - p ^ =0 (113) 

Let {4>n} be the solutions of this differential equation satisfying the boundary condition 
finiPoo) = and A n the corresponding eigenvalues. The eigenfunctions form a complete set 
of orthogonal functions for the scalar product 

rpoo 

(fg) = / f(p)g( P )27T P dp (114) 







This system can be further normalized, i.e (</> n </> m ) = S nm . Any function f(p) satisfying the 
boundary condition (109) can be expanded on this basis, and we have the formula: 

f(p) = ( 115 ) 
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In particular: 

$(P ~ Pa) = ^2 2vrp o 0n(po)0n(p) (116) 



The general solution of the problem (|107| ) ( |109| ) can be expressed in the form 

,(p, t)=J2 A n e- x ^e-P 2 / 2 n (p) (117) 



w( 



where the coefficients A n are determined by the intial condition (108), using the expansion 



( 116 ) for the 5- function. We obtain: 

w(p, t) = e ~(P 2 -^)/ 2 e- x " r Mp)MPo) (H8) 

n 

Using the foregoing solution for w, the probability that a particle initially located at po 
leaves the vortex between r and r + dr is [see Eq. (|9l])1: 



p(po,T) = -ir Poo e (p - p ° )/2 

^ e " A "7(fe)^W ( 119 ) 

The probability Q(po,r) that the particle has left the vortex at time r is therefore [see Eq. 
©]: 

Q( /0o ,r) = -vrp ooe -^- p o)/2 
V A n dp 

Finally, to obtain Q(t), we have to take the average of the foregoing expression over the rel- 
evant range of pq. To make the solution explicit, it remains to determine the eigenfunctions 
4> n and eigenvalues \ n of the quantum oscillator. 

4.4 The ground state of the quantum oscillator 

The dependance of the eigenvalues X n on the size of the "box" can be obtained from a 
procedure developed by Sommerfeld in his studies of the Kepler problem and the problem 
of the rotator in the quantum theory with "artificial" boundary conditions (Sommerfeld & 
Welker 1938, Sommerfeld &: Hartmann 1940). This was used later on by Chandrasekhar 
(1943) to determine the rate of escape of stars from globular clusters, from which our study 
is inspired. 

With the change of variables 

x = p 2 (121) 



d 2 6 dd> I 1 A x\ 



equation (|113j) becomes: 
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Another change of variables: 

cf> = fe-*l 2 (123) 
transforms equation (122) into Kummer's equation 

dx z ax 2 

In the case of a "free" oscillator (p^ — * oo), it is well-known that the eigenvalues are 
X n = 2n (with n = 0, 1, ...) and the eigenfunctions are proportional to Laguerre polynomials 
/„ = ALn ■ In a bounded domain (p^ < oo), the eigenvalues will be different but if poo is 
sufficiently large, we expect that the difference will be small. Therefore, we expect A n ~ In. 
Accordingly, for values of r of the order of unity, or greater, the first term in the series 
(120) will provide ample accuracy. Thus 

Q(Po,t) ~ -vrpooe-^--^/ 2 

(l-e~ x n^r(Poo)MPo) (125) 
A dp 

We have therefore to determine the ground state of our artificially limited quantum oscil- 
lator. To that purpose, we expand / in a series: 



f = Y,a n x n (126) 



and substitute this expansion into equation Q124| )- Identifying term by term, we obtain the 
recursion formula: 

n- * 

an+1 = ( n + l)2 Q » ( 12? ) 

which is easily reduced to 

1 2/ -«o (128) 



(n _l_ |)(n-2 - A).,.(-|) 
(ra!) 2 



We know already that the eigenvalue Ao of our artificial quantum oscillator is very small. 
Keeping only terms of order Aq in the products (128), we obtain: 



(n _l)( n _2)...l(-lf) 



(n!) 



P|2 



2 a 



(n!) 2 2 u n\n 2 
The ground state function can therefore be written 



{n- 1)! A a A 

«o = j — 7T (129) 



fo(x) = a (l - yxW) (130) 



with 

xfc) = A- = Ei(x) -Ins -7 (131) 



-oo 



, n\n 

n=l 
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where ^ 

Ei(x) = P [ e -dt (132) 

is the Exponential integral and 7 = 0.577215... the Euler constant. The function x{ x ) 



is plotted on figure 11. The fundamental eigenvalue Ao is determined by the condition 
fo(plc) = 0, i.e 

^0 = / o \ (133) 

x(PSo) 

Returning to the original eigenfunction <fo (p) , we have established that 

MP) = ao(l - yX(p 2 ))e V/2 (134) 
where ao is a normalizing factor. The final result can therefore be expressed in the form: 

Q(r) = A(l - e- AoT ) (135) 

with 

A = -^e-^^( Poo )e^M P0 ) (136) 
Ao dp 

In the expression for the amplitude, the bar indicates an average over the relevant range of 
Pq. In practice, A is close to unity but its precise value determines to which accuracy the 
approximation consisting in keeping only the dominant term in the series ( |120| ) is justified. 
For sufficiently large poo's, this will always be the case, so we shall take 

Q(t) ~ 1 - e~ Aor (137) 

This expression is consistent with the physical condition Q(t) — > 1 as r — > 00. 

4.5 Application to the solar nebula 

Let Mo denote the total mass of particles present in the vortex at time t = 0. If we assume 
no renewal from the outside then, on account of turbulent fluctuations, some particles will 
escape the vortex and, according to equation (137), the mass will decay like: 



M(t) = M e~ Xot/tcapt (138) 
The "evaporation" will take place on a typical time 

tesc — T tcapt (139) 

Ao 



where Ao is given by equation ( 133 ). Note that both Ao and t capt depend on the friction 
parameter £ of the particles. For very light or very heavy particles, p^ — ► and Ao — 2/p^. 
With the definition ( gOp we find: 

tesc = ^ (140) 
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Figure 11: The function x( x ) 



in good agreement with the result ( |101[ ) obtained when the drift term is ignored. This shows 
that equation (|139| ) can be used with good accuracy for all types of particles. Substituting 
equations (]73|)(|74|) and (|33|) (|34|) into equation (|139|), we obtain explicitely: 



11 /VtTT^ 
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(light particles) 



(heavy particles) 



(141) 



(142) 



16a V O 

where \ is the function defined by equation (|131| ). When £ — > oo, equation (|141| ) tends to 



fit. 



4a 



200 



(143) 



Very light particles are not concentrated in the vortices and they diffuse away after only 
~ 30 rotation periods. This is the minimum escape time. By contrast, for optimal particles 
with £ ~ Q, the concentration reduces considerably the evaporation and we find 



4^ / 3y^+T 

3n X \ 16a 
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61 



(144) 



The particles are so much concentrated in the vortices (Id ~ O.li?) that the turbulent 
fluctuations are not sufficient to allow them to reach the edge of the vortex. Therefore, 
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Figure 12: Escape time of the particles as a function of their friction parameter. 



their escape is completely negligible. Finally, for £ — > 0, we have 



4a 



(145) 



For very heavy particles, the escape time goes to infinity because the particles are not 
coupled to the gas and therefore not affected by diffusion. Recall, however, that our study is 
not well suited for very heavy particles. As discussed in section |2.3| , these particles are more 
likely to cross the vortices without being captured. The curve Qt esc versus f/fi is potted 
on figure [l^ (see also figure 13 for the dependance of flt esc on the size of the particles). 
The asymptotic regimes (143) and ( |145| ) are obtained for particles with £/f2 > 153 and 
< 4 IP -5 respectively. For these particles, Id > R (see inequalities fl7q) and (|7q)), so 
there is no concentration. Therefore, the asymptotic limits ( |143j ) and ( |145| ) agree with the 
results 



( |10 2D and (1P3) obtained when the drift term is ignored. 
These results assume that there is no renewal of the particles in the vortices. If we now 
take into account a capture process like in section |2.4j , we are led to consider the balance 
equation: 



dM 



1 



t. 



M 



(146) 



The first term in the r.h.s accounts for a flux of particles inside the vortex (see equation (|37|)) 
and the second term for an exponential decay of the particle number due to "evaporation" 
(see equation (|138|)). 
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Figure 13: Escape time of the particles as a function of their size at r 
(p s = 2g/cm 3 ). 



IA.U and r = 5A.U 



For particles which can experience gravitational collapse (see inequalities (^)(^3|) and 
I)), the escape time is much longer than the vortex lifetime (see figure |i~2| ). In that case, 
evaporation can be neglected and equation (146) reduces to equation (37). The maximum 
mass captured by the vortex is determined by its lifetime and the results of sections 2.5 



and |T4| are unchanged. 

Consider, however, the idealistic situation when tnf e — > oo (recall that we don't know 
precisely the value of tnf e ). For sufficienlty large times (t > t esc ), the vortex will achieve a 
stationary distribution of dust particles obtained by setting dM/dt = in equation ( |146| ) . 
This gives a maximum mass 



-M nt esc)R 2 f(0 



(147) 



which is now limited by the evaporation time (compare equation (147) with equation (38)) 
The density enhancement in the vortices (taking into account the concentration effect) is 



a 



vort 



(148) 



as represented on figure |l4[ We find that even if the vortex had an infinite lifetime, particles 
with friction parameter £/Q > 30 cannot trigger the gravitational instability. This result 
implies (see figure ||) that subcentimetric particles cannot form planetesimals even in the 
most optimistic case. Sticking processes are necessary to produce larger particles. 
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Figure 14: Enhancement of the dust surface density in the vortices as a function of the 
friction parameter when evaporation due to small scale turbulence is taken into account 

(tlif e -> OO). 
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5 Conclusion 



This article follows the works of Barge & Sommeria (1995), Tanga et al. (1996), Bracco 
et al. (1999) and Godon & Livio (1999c) concerning the interaction between dust particles 
and large-scale vortices in the solar nebula. The originality of our study was to start 
from an exact and realistic solution of the fluid equations and to provide analytical results 
and relevant parameters for the trapping process. Moreover, we have considered for the 
first time the effect of small-scale turbulent fluctuations on the motion of the particles 
and determined an explicit expression for the escape time by solving a problem of quantum 
mechanics. These theoretical results have been formulated in the context of Keplerian disks 
and planet formation but they can clearly have applications in other fields of astrophysics 
and geophysics, for example the transport of polluants in the Earth's atmosphere. 

The vortex scenario provides an attractive mechanism to form planetesimals on very 
short time scales. This does not contradict the overall picture of planet formation which 
has been developed over the last decades. On the contrary, it fills the gap between two 
domains which were difficult to connect: sticking processes are still necessary to produce 
centimetric particles and collision between planetesimals is of course the main engine of 
planet's growth. In between, the vortex scenario should come at work to facilitate, at some 
prefered locations of the disk, the Safronov-Goldreich-Ward instability which was proposed 
initially for forming planetesimals. 

The predictions of the vortex model are remarkably consistent with the structure of 
the solar system. The capture process is optimal at two prefered locations of the nebula 
which correspond, for relevant sizes of the particles (i.e some centimeters), to the position 
of telluric and giant planets. The transition between the two groups of planets happens 
to coincide with the passage from the Stokes to the Epstein regime where the gas drag 
law changes. The asymmetry between the two optima may be ascribed to the size of the 
vortices which are bigger in the outer zone. 

Of course, the results discussed here rely upon the existence of vortices in the disk. Their 
presence in the solar nebula is reasonable due to the ubiquitous appearance of vortices in 
rotating flows and two-dimensional turbulence. However, numerical simulations and even 
laboratory experiments are necessary to ascertain their existence in Keplerian flows. A 
first step was undertaken by Bracco et al. (1998,1999) and Godon & Livio (1999a,b,c) who 
observed the formation of anticyclones developing from an energetic random initial vorticity 
field in a Keplerian flow. More work remains to be done to understand the generation 
(and maintenance) of vorticity in such disks (convection, baroclinic instability...) and to 
determine the effect of three-dimensionality (Hodgson & Brandenburg 1998) on the stability 
of theses vortices. 
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A Elliptic vortex in a uniform shear 



In his monograph, Saffman (1992) reports the existence of an exact solution of the incom- 
pressible two-dimensional Euler equation consisting of an elliptic patch of uniform vorticity 
ijj embedded in a simple shear k. Since this vortex solution is the starting point of our 
study, we establish in this appendix the condition (|l^) for its existence and construct the 
corresponding streamfunction. These results will be useful for subsequent studies. 

First, we introduce a cartesian system of coordinates and write the shear under the 
form 

u x = ny (149) 

u y = (150) 

where k is a constant (k = |f2 for the Keplerian shear considered in section ^^)- The 
associated vorticity is u ex t = — «. 

Inside the vortex, the velocity field can be written 

q 2 

u x = -— — ^uy (151) 
l + q z 

' T ux (152) 



where q = a/b is the aspect ratio of the elliptic patch (a and b denote the semi-axis in the 
x and y directions respectively) and u the vorticity. We can check that the fluid particles 
move at constant angular velocity j^-z^ along concentric ellipses with aspect ratio q. 

For an incompressible two-dimensional flow, it is convenient to introduce a streamfunc- 
tion if) defined by u = — z A Vip (where z is a unit vector normal to the flow). For a given 
vorticity field, the stramfunction is solution of the Poisson equation 

AV> = -u (153) 

Inside the vortex, we find that 

^ = -\y^(* 2 + fy 2 ) (154) 

where we have taken, by convention, ifj = at the centrer of the vortex. On the vortex 
boundary, the streamfunction is constant with value: 

Outside the vortex, we have to solve the Poisson equation 

Aip = k (156) 

with boundary conditions at infinity 
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These boundary conditions insure a continuous matching with the shear ( 149[ ) fll5C)| ) at large 
distances. Introducing the decomposition 

^ = cf>+~Ky 2 (158) 

the problem ( |156[) (|15?1 ) is equivalent to solving the Laplace equation 

A(f) = (159) 

with boundary conditions at infinity 

^ -> 0, ^ for x, y oo (160) 
ax ay 

At this stage, we find it convenient to use elliptic coordinates 

x = c cosh £ cos 7/ (161) 

y = c sinh £ sin 77 (162) 

with < £ < oo and < r/ < 2n. The vortex boundary is an ellipse with parameter £o 
satisfying 

jl + o ^27 = 1 (163) 



c 2 cosh 2 £o c 2 sinh 2 £o 
This relation determines the semi-axis a and b of the ellipse in terms of £o and c: 

a = ccosh£o 6 = csinh£o (164) 

Alternatively we have 

tanh£ = - = - c 2 = a 2 -b 2 (165) 
a q 



In terms of elliptic coordinates the Laplace equation (|159| ) has the simple form 

^2+^2=0 (166) 



This equation, with the boundary condition ( |160| ), is solved easily and we obtain outside 
the vortex 

j +oo 

"0 = 2 K V 2 + B Z + Y1 An€r< cos( - nr l + 7») ( 16? ) 

n=0 

where -B, ^4 n and 7 n are some constants. At large distances, £ — > lnr and r/ ^ 9 (where r 
and are polar coordinates) and we recover the usual multipolar expansion of the stream- 
function. The condition that the vortex boundary is a streamline destroys most of the 
terms in the series ( |167D . There remains only 

tj) = — ^— sinh 2 £(1 — cos(2r/)) 
Kb 2 

+B£ + A + — e 2(f °- ?) cos(2r/) (168) 
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In addition, the continuity of tp on the vortex boundary requires that Aq and B be related 
by 

- l-^a 2 = ]Kb 2 + B^ + A (169) 
2 1 + q z 4 

We must also satisfy the continuity of the tangential velocity ^ at the contact with the 
vortex. To that purpose, we need first to express the streamfunction ( 154 ) inside the vortex 
in terms of elliptic coordinates. We find: 

ip = — — 2-c 2 (cosh 2 £ cos 2 i] + q 2 sinh 2 £ sin 2 rj) (170) 

2 1 -\~ q 

Then, after some algebra, we find that the continuity of the velocity is satisfied provided 
that 

uo 1 + q z 

and 

B = -^(k + u)qb 2 (172) 

The relations (172) and (169) just determine the constants Aq and B appearing in the 
expression ( |16Sj ) of the streamfunction. In addition, the condition (|171| ) must be satisfied 
for a solution to exist. In a given shear, equation ( |171| ) imposes a relation between the 
vorticity and the aspect ratio of the vortex. 

Regrouping all the results, we find that the streamfunction can be expressed inside the 
vortex (£ < £ ) by: 

Kb 2 

ipin = (q + l)(cosh 2 £ cos 2 r] + q 2 sinh 2 £ sin 2 r\) (173) 

2q 

and outside the vortex (£ > £o) by: 

VW = - 1) sinh 2 £(1 - cos(2r ? )) 

L 1 ^ 9 + 1 , 1,2 9+1/ 1 C N 

4 g — 1 2 g — 1 

+ ^ e 2«o-O cos ( 2r? ) (174) 

B Gravitational instability of a turbulent rotating disk 

In this appendix, we derive an instability criterion for a turbulent rotating disk in the 
approximation where the disk is a sheet of zero thickness with constant mean density a 
and constant angular velocity 17. Our study is inspired by the work of Chandrasekhar (1951) 
concerning the stability of an infinite homogeneous turbulent medium. To our knowledge, 
the results presented in this appendix are new. 

The analysis is easiest if we work in a rotating frame of reference. The equations of the 
problem are provided by the equation of continuity, the equation of motion and Poisson 
equation: 
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dt dxj J ' Q Xi 



+cr0 2 Xi - a^- (176) 
ox 



A$ = AitGa5{z) (177) 

The symbols have their usual meaning and uj_ is the vector u rotated by \. The velocity 
field Ui = (u x ,u y ) is purely two-dimensional and, in equations (|175[)(|176| ), there is sum- 
mation over repeated indices. In equation (177), the 5-function insures that the disk is 
infinitely thin. 

We assume that the disk is turbulent and write 

a = a + 5a, p = p + 5p, $ = ¥ + 5$ (178) 

where a, p and <3? are certain constants. With the further assumption that the disk is 
barotropic, i.e p = p(a), we have 

dp 9 da ,„ 

SH&J (179) 

where c s = (dp/ da) 1 ' 2 denotes the velocity of sound. Then, we invoke Jeans swindle (see, 
e.g, Binney & Tremaine, 1987) to eliminate the centrifugal term, i.e we write 

— - Q 2 x t = -— 180 

OXi OXi 

This process is permissible if we assume that the centrifugal force is balanced by a grav- 
itational force that is produced by some unspecified mass distribution (recall that for an 
infinite uniform disk, V<£ is necessarily vertical and cannot, by itself, compensate the cen- 
trifugal force). In our situation, the external force is provided by the sun's gravity. With 
the further approximation 

85$ . ^ 85$ 85$ , . 

a—- = [a + 5a — ^a— 181 

8xi 8xi 8xi 

the equations of the problem become 

£ + £(„,)_<> (182) 

d . . 8 o8a _85$ . 

di i<7Ul) + = ~ Cs 8x~ ~ 2naU±l ~ a ^x~ (183) 

A5$ = ^G5a5{z) (184) 

From the continuity equation ( |182|) , we readily establish that 

P^ = ~i^-^^ (185) 
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where a and a' are the values of the surface density in x and x' respectively. We introduce 
the correlation function of the density fluctuations, defined by: 



C(f , t) = 5a5a' = (a -a) (a' -a) = an' - a 1 (186) 

In homogeneous and isotropic turbulence, C(£, t) is a scalar function depending, apart from 
time, only on the relative distance £ = |x' — x| between the points under consideration. 
Similarly, we define the quantity 



Li(£, t) = cr'aui = -aa'u^ (187) 

and set 

f = x' - x (188) 
Remembering that d/dx^ = —d/dxi = d/d£i, we can rewrite equation ( [185[) under the form 

8(1 

— = 2V € L (189) 

An equation of motion for L(£,i) can be derived in the following manner. Prom equations 



( 182 ) and (|183| ), it is easy to establish that 



d —. d — ; r d — 

—a'aui + —-aa'uiu': + - — a'auiUj 
at ox- J oxj 



2 da'a da' 5^ . . 

-c^ a— 190 

OX; OX; 



or, with more convenient notations 



9Li d — r 2 9C _dtp 

— + Wj aa'(u iUj - u t u 3 ) = c s — + a— - 2QL ±l (191) 



where we have written 



if; = a >5<S> = 5a>5<S> (192) 
The correlation function ip is related to C by the Poisson equation 

AV> = ^GC5{z) (193) 

Equations ( pl9|) ( |l9l| ) and (p^ ) are perfectly general but the system is not closed since 
the evolution of L(£, t) involves the fourth-order correlation functions aa'uiu'j and aa'uiUj 
which are not known. To go further, we shall suppose that these fourth-order correlations 
can be expressed in terms of the second-order correlations as in a joint gaussian distribution. 
For our purposes, this approximation is reasonable and should provide ample accuracy. 
Then, by virtue of Wick's theorem, we have: 



liUj = oo' UiUj + aui a'Uj + cruj a'u 



<Ja'^v?5 l] (194) 
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a a' UiUj + aui a'u'j + au'j a'ui 



cru'i au'j 



(195) 



With the foregoing substitution, equation (191) becomes 



dL 
~di 



0_, r 



u 2 dC 

dib „ 
- 2QL ±i 



(196) 



For sufficiently large values of £, the terms in the velocity correlations in equation (|19 
will become negligible and equation (196) will tend to 



dL 

~dt 



( c s + y )V ? C + aV^ - 20 A L 



(197) 



Equations ( p.89| ) ( |197| ) and ( |193| ) are mathematically similar to the linearized equations which 
appear in the usual problem of Jeans instability for a thin rotating disk (see, e.g, Binney 
Sz Tremaine 1987). However, their physical meaning is different since equations ( |189| )(197) 
and ( 193 ) have been derived explicitely from the analysis of the correlations in a turbulent 
medium. 

These equations admit sound waves of the form 



C 
L 

ip = ipe 



iujt J (hOe- klzl 



(198) 

(199) 
(200) 



where Jo is Bessel function of order zero. The solution ( |200| ) satisfies the Laplace equation 
£±ip = for z / [see equation (|193|)]. In the plane z = 0, we have ip = tpe~ llvt Jo(k( t ). 



Equation (193) implies that ip must be related to C in a special manner. To see that, 
we integrate the Poisson equation from — £ to (where £ is a positive constant) and let 
C - 0. Since g$ and ^ 



d r 2 



are continuous at z = 0, but -wX is not, we have 



lim 



< d 2 4> 



c 

4ttGC 



dz 2 



dz 



lim 



5{z)dz 



dip 
dz 

: 4ttGC 



C 

-C 



Hence —2kif) = 4irGC, or 



2t:GC 
k 



Substituting for C, L and ip in equations ( |189| )( ^97| ) and (193), we obtain 

1 d 



iujC 



2-^-(£if) 



(201) 



(202) 



(203) 
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/ 2 u 2 .dC _dip 



where £, r\ are polar coordinates. Eliminating L, n between ( |204| ) and fl205|) yields 

■ 2 4ft 2 . , „ ^.dC 



Substituting the foregoing expression for into equation (|203| ), we arrive at 

2 

(4ft 2 - uj 2 )C = 2(c 2 + — )A C C + 2ctA^ 



(204) 
(205) 

(206) 



(207) 



Using the expressions (198) ( pOOj ) ( |202j ) for C and -0 an d remembering that Jo(k^) is an 
eigenfunction of the two-dimensional Laplacian operator with eigenvalue — k 2 , we obtain 
the dispersion relation 



An 2 + 2 



k 2 (c 2 s + ^—)- 2-KGak 



(208) 



This result differs from the usual dispersion relation (77) in the occurence of c 2 + \ in 
place of c 2 and by a factor 2. Similar differences with the usual Jeans dispersion relation 
were noticed by Chandrasekhar (1951) in his analysis of an infinite homogeneous turbulent 

medium. 

The function uj 2 (k) is quadratic in k with a minimum at k c = irGa/(c 2 + ^-). The disk 
will be unstable to some wavelengh k~ l provided that u> 2 (k c ) < 0, i.e: 




2 cj + 



< 



(209) 



This is the generalisation of the Toomre instability criterion fl7q ) in the case where the disk 

2 2 

is turbulent. For a real disk, with finite thickness H, ^- should be replaced by ^ since 
small-scale fluctuations are basically three-dimensional. When the turbulent dispersion 
c turb = V^ii 2 dominates over the sound velocity, as it is the case in our study, we obtain 



Cturb ^ 



n 



(210) 



This result justifies the procedure used in section ^4| of replacing the sound velocity occuring 
in equation (|7^) by the more relevant turbulent dispersion of the particles. However, in more 
general situations, the turbulent dispersion needs not be large compared to the velocity of 
sound and the criterion (|209|) should be used. 
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